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Electric field representation of pulsar intensity spectra 
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ABSTRACT 

Pulsar dynamic spectra exhibit high visibility fringes arising from interference between scat- 
tered radio waves. These fringes may be random or highly ordered patterns, depending on the 
nature of the scattering or refraction. Here we consider the possibility of decomposing pulsar 
dynamic spectra - which are intensity measurements - into their constituent scattered waves, 
i.e. electric field components. We describe an iterative method of achieving this decompo- 
sition and show how the algorithm performs on data from the pulsar B0834H-06. The match 
between model and observations is good, although not formally acceptable as a representa- 
tion of the data. Scattered wave components derived in this way are immediately useful for 
qualitative insights into the scattering geometry. With some further development this approach 
can be put to a variety of uses, including: imaging the scattering and refracting structures in 
the interstellar medium; interstellar interferometric imaging of pulsars at very high angular 
resolution; and mitigating pulse arrival time fluctuations due to interstellar scattering. 

Key words: pulsai-s: general — ISM: general — scattering — turbulence — techniques: 
interferometric 



1 INTRODUCTION 

In the last few years there has been much progress made in under- 
standing the interference fringes which are manifest in pulsar dy- 
namic spectra. Representing the data as power spectra ("secondary 
spectra") of the dynamic spectra demonstrated the underlying sim- 
plicity of the complex organised fringe patterns which are some- 
times seen: commonly it is found that the power is concentrated on 
parabolic loci in the transform domain (Stinebring et al 2001). This 
property is now understood in the following terms (Stinebring et al 
2001; Walker et al 2004; Cordes et al 2005). The dynamic spec- 
trum is the electric field intensity, I(i',t), as a function of radio- 
frequency, u, and time, t; the conjugate (Fourier) variables are the 
delay, r, and Doppler-shift, uo, of the scattered waves. The purely 
geometric component of the delay is proportional to the square of 
the scattering angle, whereas the Doppler-shift is proportional to 
one component (the component parallel to the effective transverse 
velocity vector) of the scattering angle, leading to parabolic rela- 
tionships between r and u. 

The well-defined parabolae which are often seen further re- 
quire that the scattering material not be distributed along the line- 
of-sight, but instead must be concentrated into a thin "screen" 
(Stinebring et al 2001). In some cases the data also suggest that 
the scattering is highly anisotropic (Walker et al 2004; Cordes et al 
2005). 

To date the interpretation of observed dynamic/secondary 
spectra has proceeded by forward theoretical modelling; in other 
words models are constructed and the predicted secondary spectra 



are compared with the data. This approach has yielded important 
insights, but it is fundamentally limited in that a detailed match 
to the data is not practicable — only the general characteristics 
can be considered. This is a severe limitation because dynamic 
spectra may contain a great wealth of detailed information in the 
~ 10® independent frequency-time measurements which can be 
routinely recorded. Consequently we are motivated to attempt in- 
verse modelling of dynamic spectra, whereby the electric fields are 
deduced by modelling the observed intensity distribution. If the 
total electric field is U{v,t), then I{iy,t) — U''{u,t) U{y,t), so 
\U{v,t)\ = ^ I{u, t) tells us the field amplitudes, but the phases 
remain a priori unknown. This is an example of a "phase retrieval 
problem"; such problems are common in the optics literature (e.g. 
Fienup 1982; Elser 2003; McBride, O'Leary and Allen 2004), al- 
though we were unaware of that resource until the present work 
was largely complete. 



This paper presents an iterative approach to solving the inver- 
sion problem; it is not the only approach, nor is it necessarily the 
best approach, but it has been at least partially successful. The out- 
line of the paper is as follows: in the next section we describe the 
physical model we have adopted and derive the mathematical ba- 
sis for our inversion algorithm; our attempt at implementing this 
algorithm is described in §3, while in §4 we show how this im- 
plementation performs on real data; possible applications of these 
techniques are discussed in §5. 
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2 FOUNDATION OF THE INVERSION 

As radio-waves propagate through the Galaxy they are subject to re- 
fraction and scattering by inhomogeneities in the ionised interstel- 
lar medium (e.g. Rickett 1991; Narayan 1992). If the radio-waves 
originate from a pulsar then the size of the source is small, and the 
coherent patch correspondingly large, so that essentially all pairs of 
scattered/refracted waves yield high-visibility interference fringes, 
regardless of the spatial separation of the scattering centres. Con- 
sequently the combination of scattered waves can be well approxi- 
mated simply by addition of the various electric field components. 
(We consider departures from this point-like source approximation 
in §5.3.) Denoting the total electric field as U (v, t), and each of the 
discrete scattered waves as Uj {v, t), we have 

U {v, t) = ^ Uj {v, t) = ^ Uj exp [2ni{i'Tj + u)jt)] . (1) 

Equivalently, we can write this relationship in the Fourier Trans- 
form domain, {v, t) — » (r, w), as 

U (r, = ^ 5j S{t - Tj) 5{oj - u)j), (2) 

where S denotes the Dirac delta-function, and the Uj are complex 
constants. The intensity observed at radio frequency v (relative to 
band centre) and time t (relative to the central epoch of the obser- 
vation) is I(v,t) = U*{v,t) U{i^,t). Given an observed I{u,t), 
we want to determine U {v, t). More precisely, because the conju- 
gate variables (r, lS) have clear physical interpretations, as the de- 
lay and Doppler-shift of the scattered waves, we wish to determine 
U(t, uj). In this model the individual scattered waves which make 
up U are completely specified by their delay (jj), Doppler-shift 
{uij), amplitude and phase (the single complex number Uj\ 

We have taken an iterative approach to solving for U , given 
/, as we now describe. Given a model electric field, C/q, we can 
compute the corresponding intensity pattern. In general this will not 
match the data exactly and we want to improve the model. Suppose 
our model differs from the true electric field by an amount 5Uo, we 
can then write 

I{v,t) = U*U = US Uo + 5US Uo + UoSUo + SUS 5Uo . (3) 

(Where there is no ambiguity we will henceforth not make explicit 
the independent variables.) Introducing the residual between model 
and data, TZoii^, t) := I — UqUo, we arrive at 

7^o = 5U; Uo + Uo 5Uo+5U55Uo. (4) 

Providing the existing model is a good one <S we can 

neglect the last term in this expression and estimate the quantity 
SUo from the resulting linear approximation. 

Now suppose that 5Uo is dominated by a single scattered 
wave, so that our task reduces to a determination of the properties 
of that wave, then we can make the approximation 

5Uo {v, t) = 5Uo exp [2m{vT + u>t)] , (5) 

for some particular (but unknown) values of (r, to), where 5Uo is a 
complex constant. Dropping the final term in equation 4, multiply- 
ing by Uoiiy, t) exp [— 27ri(z/r -|- ujt)] and integrating yields 

J dt df7?.o?7o exp [—2m{i'T + ujt)] = 

5Uo J dtdf Uq exp [—A-niiyT + ujt)] 

+ 5Uo UtAvUoUo- (6) 



Although the first two terms on the right-hand side of equation 4 
are comparable in magnitude, their counterparts in equation 6 are 
not; generally the term in 5Uo is expected to dominate the term in 
5Uq , because only a small subset of the power in Uo contributes 
to the coefficient in SU*), whereas (by Parseval's theorem) all of 
the power in {/o contributes to the coefficient of dUo. Noting that 
JdtdvUoUo ^ /dtdi^/(z/,t), because |5C/o I <C |i7o|, we see 
that we can choose J dt df Uq ?7o = 1 at the outset, by appropri- 
ately normalising the dynamic spectnmi. Neglecting the conjugate 
image (i.e. the term in SUq) then leads to the simple form 

SUo = nolJo. (7) 

This result gives us an estimate for the (complex) amplitude of the 
scattered wave which is missing from the model, as a function of 
the assumed delay and Doppler-shift of that wave. This result has 
been derived under the assumption that MJq is dominated by a sin- 
gle scattered wave, so that the appropriate choices of delay and 
Doppler-shift are those values for which the modulus of the right- 
hand-side attains its largest value. 

The result just derived shows how, given a model electric field 
[/(), we can improve that model by adding a single scattered wave 
component. Moreover the process can clearly be iterated, so that 
the restriction embodied in equation 5 - i.e. that 5Uo can be ap- 
proximated by a single plane wave - is unimportant: the full spec- 
trum of scattered waves can be built up iteratively. This is done by 
adding one new component to the reference model, Un, on each it- 
eration. In its overall structure this procedure is similar to that of the 
CLEAN algorithm (Hogbom 1974) which is commonly employed 
in aperture synthesis imaging in radio astronomy; that algorithm 
inspired the one presented here. 

We need to assume a model as a starting point for an iterative 
solution; in the absence of any specific a priori information about 
the scattered waves, the most sensible choice is a zero-frequency 
(r, = 0) plane wave model. This choice may correspond, phys- 
ically, to an unscattered wave. However, we note that the interfer- 
ence fringe properties depend on delay/Doppler differences, so our 
model_predicts the same dynjtmic spectrum for the scattered wave 
field Ui{t, uj) as it does for U2{r, uj) = Ui{t + to,oj + ujo). This 
point can be recognised most easily if we write down one of the 
interference terms contributing to / = U*U: 

uluj cxp[2Tvi{iy{rj — Tk) + iujj — ujk)t}]. (8) 

This degeneracy means that there is no information in the dy- 
namic spectrum on the absolute delay and Doppler shift of the 
scattered waves, just as there is no information on absolute phase 
(UfeUj = |Mfe I |wj I exp[i((;/>j — cjjk)]), and the origin in our model U 
is arbitrary. 

2.1 Weak scattering limit 

To our knowledge only one previous attempt has been made to de- 
rive the scattered wave spectrum from a recorded dynamic spec- 
trum: B.J. Rickett (personal communication 2003) derived a solu- 
tion {uj ,Tj,u)j}'m the weak scattering limit. In this limit all of the 
scattered wave components are of very low amplitude in compari- 
son with the unscattered component, and to obtain a solution all of 
the cross-terms u*uk with j, fc 7^ are neglected. This approxima- 
tion has a simple correspondence with the approach we have de- 
scribed. Because the scattered components are all very weak, there 
is no need to keep refining the model electric field, Uo, and conse- 
quently it is not necessary to iterate the solution: we simply take all 



of the field components SUo returned by equation 7 for the start- 
ing model, Ug = l.JThis yields the solution in the weak scattering 
limit: ?7 = ?7o + 6Uo = I{t,ui). Our approach, by contrast, em- 
ploys an unscattered reference wave {Uo = 1) only for the first 
cycle of the iteration process; all subsequent cycles differ from the 
weak scattering limit by including the identified scattered compo- 
nents in the reference model, and in general f/o 1 in equation 7. 



3 AN IMPLEMENTATION OF THE ALGORITHM 

To evaluate the procedure described in the previous section we first 
tested it on a simple synthetic dynamic spectrum made up of four 
scattered wave components plus a small amount of pseudo-random 
"noise". Whereas the derivation in §2 makes no particular assump- 
tion about the sampling of I{i', t), our synthetic dynamic spectrum, 
and the real data discussed in §4, were sampled on a regular grid 
in radio-frequency and time. In Fourier space, our representation 
of these data need only employ wave components with conjugate 
variables (delay and Doppler-shift) which lie on a uniformly spaced 
grid satisfying the Nyquist sampling criterion; our algorithm em- 
ploys such a grid. 

We found that the algorithm did indeed locate (in delay- 
Doppler space) the electric-field components which we knew to 
be present, but it did not make accurate estimates of the (complex) 
field amplitudes when the components were first identified. This 
is not surprising because the procedure described in §2 employs 
simplifying assumptions, so equation 7 is not expected to give a 
very accurate estimate of the component amplitudes. This property 
would be unimportant if subsequent iterations refined the values of 
the component amplitudes so that the algorithm converged on the 
correct solution. In fact the procedure generated spurious compo- 
nents, instead of converging on the correct values of the actual field 
components, thus smearing power around in delay-Doppler space. 
This behaviour yields results of low dynamic range and is unac- 
ceptable. 

Because of the poor performance of the procedure described 
in §2, we decided to refine the algorithm slightly. Instead of simply 
accepting the largest Fourier coefficient in equation 7 as the appro- 
priate estimate of the wave amplitude, we refined the estimate by 
adjusting the wave amplitude and phase so as to achieve a least- 
squares fit of the revised model to the dynamic spectrum. In this 
scheme the estimate provided by equation 7 is used to fix t and 
uj for the new component, as before, but the (complex) amplitude 
indicated by equation 7 is used only as the starting point for the 
least-squares minimisation. By fitting to the data we sidestep many 
of the concerns which might otherwise arise in cormection with 
the approximations made in deriving equation 7 — linearisation, 
the single component approximation, and neglect of the conjugate 
image. In short the revised procedure should work well providing 
only that the largest Fourier component of TZqUo is also the largest 
Fourier component of the difference between the reference field and 
the actual field {SUa = U — Uo). We found that this algorithm did 
indeed perform well, yielding reconstructions which matched the 
input as closely as possible (i.e. to within the "noise"), and did not 
generate any spurious components beyond what would be expected 
given the "noise" level of the "data" — i.e. the algorithm did not 
limit the dynamic range of the results. In summary, the successful 
algorithm follows the procedure specified in §2, but at each itera- 
tion the amplitude and phase of the new electric field component 
(identified via equation 7) is determined by least squares fitting to 
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the data. The algorithm is described in detail in the following sec- 
tion. 



3.1 Details of the algorithm 

The algorithm was implemented in IDL, a high-level, commercial 

software package which provides straightforward array manipula- 
tion and data display. The main elements of the algorithm are: 

(i) The data are preconditioned; this involves three steps. First 
we remove the spectral profile, B{u), which is imposed by the 
bandpass filter; this is part of our routine data reduction process, 
with -B(i') determined from a calibration data set. Secondly, the ef- 
fects of intrinsic fluctuations in the pulsar's flux are removed, in so 
far as possible; these fluctuations manifest themselves as a purely 
temporal modulation, f{t), with a white noise spectrum. We dealt 
with this modulation by forming the frequency-averaged intensity 
at each timestep; we then filtered out the high-frequency compo- 
nents in this quantity, by mutiplying by a Gaussian function of full- 
width-half-maximum equal to one half of the Nyquist frequency, 
and transformed back to yield a smoothed version of the average 
intensity; finally we made our estimate of f{t) from the ratio of 
these two quantities. Correcting the data for the effects of pulse- 
to-pulse variations is then simply a matter of dividing I{i',t) by 
f{t). It must be acknowledged that this prescription is not ideal, as 
it does not completely eliminate the intrinsic flux variations, and 
it also attenuates the wave interference structure slightly. The final 
step is to normalise to unit mean intensity. 

(ii) The starting model for Uo is taken to be a wave of unit 
amplitude and zero phase, with t, w = 0. 

(iii) A new scattered wave component is then identified as the 
largest component in equation 7 which has r ^ 0. Physically we 
expect that all delays should be non-negative relative to the unscat- 
tered wave; the algorithm can in principle differentiate between 
waves of positive and negative delay, so it should be able to dis- 
cover this property. However, if the unscattered wave is very strong, 
so that U ~ (5(t) <5(a;), then Uo ~ Uq and it is difficult for any 
practical scheme to avoid confusion between a scattered wave and 
its complex conjugate. As the latter are aU equivalent to scattered 
waves which have t < 0, there is potential for choosing spurious 
components, and in practice we found that this did indeed happen. 
To avoid this problem we imposed the restriction t > on all 
components. 

Our input data consist of measurements of field intensity on 
a regular grid in radio-frequency and time - the dynamic spectnun 
- so for these data the Fourier plane (delay-Doppler) representa- 
tion requires only components on a regular, Nyquist-sampled grid. 
Consequently we chose our scattered wave components to lie on a 
Nyquist-sampled grid in (r, cu). 

(iv) The amplitude and phase of the new component are ad- 
justed by minimising the sum of squares of the difference between 
the model and the data. Minimisation of several component am- 
plitudes/phases is undertaken simultaneously, using the "Amoeba" 
algorithm described in Numerical Recipes, which is implemented 
within the IDL software package. The 'scale' parameter used by the 
Amoeba algorithm is initially set to 0.2, when a new component is 
first introduced, and then is reduced by a factor 0.7 on each subse- 
quent iteration (i.e. every time another component wave is added to 
the model). Only components whose scale parameters are greater 
than some value (we employed 2xl0~^)are included in our least- 
squares minimisation; these criteria mean that 13 scattered wave 
components are routinely included. In addition, because the irnscat- 



4 Walker & Stinebring 



tered wave (r, tj = 0) is often very strong this component is also 
included in the optimisation, making 14 wave components in total. 

This gradual "freeze-out" of component amplitudes, as the 
Amoeba scale parameter gradually decreases, was incorporated 
into our algorithm in order that components with similar amplitudes 
could be simultaneously adjusted (optimised) in an efficient way. 
This aspect of the algorithm has an additional benefit: it helps to 
guard against the possibility of spurious parameter determinations 
arising from local, rather than global minima found by Amoeba. 
If Amoeba finds a local minimum, with correspondingly erroneous 
component amplitudes, on a given iteration, it may well find its way 
out of the local minimum to reach the global minimum on the next 
iteration cycle. Only one of the 14 component amplitudes ceases 
to be adjusted on each cycle, so the algorithm has a fair degree of 
robustness to the potential problems caused by local minima in the 
chi-squared hyper-surface. If, for some reason, the algorithm were 
to fix a component amplitude at a value which is badly in error, 
then that error may be fixed in later iteration cycles, as there is no 
barrier to putting new components in the same place as existing 
components if that is where the largest difference between model 
and data occurs. 

(v) We used the reduced chi-squared value (chi-squared per 
degree of freedom) to measure the success of our model in fitting 
the data. Each new scattered wave component that is added to the 
model should lower the reduced chi-squared, if that component is 
significant; however, once the algorithm reaches the noise level, 
the new components are (by definition) no longer significant, and 
the addition of any given component is just as likely to increase 
the reduced chi-squared as to decrease it. To reflect this change 
in behaviour we forced the algorithm to terminate when 3 succes- 
sive iterations (i.e. new scattered wave components) each caused 
the reduced chi-squared value of the fit to increase. This criterion 
causes the algorithm to terminate very quickly once the noise level 
is reached, while remaining robust to the presence of a small num- 
ber of insignificant components in the solution (i.e. the algorithm is 
not tripped up by isolated insignificant components). It is possible 
to check the significance of each wave component in the solution 
set {uj , Tj , ujj } returned when the algorithm terminates, and then 
cull insignificant components; this would be a reasonable require- 
ment to enforce, but the fraction of such spurious components is 
expected to be very small and to date we have not employed any 
culling. 



4 PERFORMANCE WITH REAL DATA 

To be useful, the algorithm must be able to cope with real data. We 
tested our code on a dynamic spectrum of the pulsar B0834-I-06, 
taken at an observing frequency of 321.0 MHz with the Arecibo 
radio telescope. We used the WAPP backend signal processors to 
record 1024 channels across atotal of 1.563 MHz of bandwidth. We 
sampled the spectrum every 4.096 ms and formed a pulse-phase- 
averaged estimate of the on-pulse and off-pulse spectra which we 
used to calculate the ON - OFF spectrum. We then averaged this 
over 10 s to form one column in the dynamic spectrum. The results 
of a 45 minute observation are shown in figure 1 in the form of 
the measured dynamic spectrum (top left) and secondary spectrum 
(power spectrum; bottom left), compared with the corresponding 
values reproduced by our model (right). From these results we can 
see that the model yields a convincing representation of the data: 
the deficiencies of the modelling are not apparent to the eye in ei- 
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Figure 1. Data (left) and model (right) for an observation of 
PSR B0834+06, in a 1.563 MHz band centred on 321.00 MHz. The data 
were taken with the Arecibo radio telescope in conjunction with the WAPP 
backend signal processing units, on MJD 53009; there are 1024 spectral 
channels, and 270 time samples, each of 10 seconds duration. The top pan- 
els show dynamic spectra, while the lower panels show the corresponding 
secondary spectra (power spectra of the dynamic spectra). Inverse grey- 
scale (black is peak intensity) is used in all cases; the transfer function is 
linear for the dynamic spectra, and logarithmic for the secondary spectra. 
The signal/noise ratio on each pixel of the observed dynamic spectrum is 
2.7, on average, yielding a dynamic range of 63 dB on the corresponding 
secondary spectrum. The model reproduces the data well. 



ther the dynamic spectrum or the secondary spectrum representa- 
tions. 

The model dynamic and secondary spectra shown in fig- 
ure 1 are of course derived quantities; the model itself is the 
set of scattered waves (electric field components) represented by 
{ujyTjyOJj}. In figure 2 we show the scattered wave amplitudes, 
|w(r, t<;)|; the roughly parabolic relationship between r and u 
is evident in this plot. These results also show how the scatter- 
ing/refracting centres are picked out very clearly in this representa- 
tion of the data as power concentrations in delay-Doppler space. It 
is beyond the scope of this paper to analyse the information present 
in the scattered wave solutions, so we will say nothing specific 
about the meaning of the results shown in figure 2; nor do we plot 
the phases of the various components. Various applications of our 
technique are, however, discussed in general terms in §5. 

In fact our scattered wave model for these data is not ac- 
ceptable, in the sense that its chi-square value is too large, given 
the number of degrees of freedom. There are 1024 x 270 inde- 
pendent pixels in the dynamic spectrum; in the model there are 
8720 scattered waves, each of which is described by four param- 
eters (amplitude, phase, delay and Doppler- shift), and the total 
over all pixels of the squared-residual between model and data is 
2.88 X 10^, in units of the noise on the dynamic spectrum. This 
corresponds to a reduced chi-squared value (Xr> chi-squared per 
degree of freedom) only slightly larger than unity: Xr = 1-19 = 
2.88 X 10V(1024 X 270 - 8720 x 4); but a statistically acceptable 
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Figure 2. The amplitudes, \uj\, of the 8720 scattered wave components 
identified by the algorithm described in §3, operating on the dynamic spec- 
trum shown in figure 1. These wave components form the basis from which 
the model dynamic and secondary spectra in figitfe 1 are derived. The axes 
are delay (t), and Doppler shift (a;). In this figure we can see that the indi- 
vidual scattered waves cluster tighdy around a parabolic locus with t oc w^. 
An inverse logarithmic grey-scale is used for this figitfe. 

model would have a reduced chi-square value much closer to unity 
(|x? — 1| ^ 0.01), so although figure 1 is impressive it is possible 
to do much better. In other words the differences between model 
and data are not due to noise alone. 

A clearer test of the performance of the algorithm comes from 
differencing the model and observed secondary spectra, as the latter 
exhibits a very high dynamic range. We find that this difference has 
a dynamic range of 47 dB, compared to the 63 dB dynamic range of 
the input data. These figures confirm that the model accounts for the 
majority of the structure in the data, but the residuals are far from 
noise-like. We expect that substantial improvements in the quality 
of the fit would be possible if the 8720 (complex) component am- 
plitudes were simultaneously optimised. Global least-squares opti- 
misation of our solution has not been attempted by us. With such 
a large number of free parameters, a global optimisation is not an 
easy task: in a linearised approach the simultaneous equations we 
are required to solve have ~ 3 x 10* non-zero coefficients. 



5 DISCUSSION 

As noted above, the algorithm we have described does not gener- 
ate a statistically acceptable description of the data shown in figure 
1. However, it was not obvious at the outset of this study that the 
type of spectral decomposition process we sought would work at 
all, so the partial success we are able to report is in fact quite en- 
couraging. We expect that much better fits to data will be achieved 
in future with further development of this technique. In part our 
confidence stems from the fact that a pulsar dynamic spectrum can 
be expected to be accurately modelled by the form I(y, t) = U*U 
(once calibration etc. is taken care of), and the job is simply to find 



U. More precisely, the job is to find the phase of U, since the am- 
plitude is known directly from I. This emphasises the importance 
of a fact mentioned in the introduction: this type of problem, i.e. 
phase-retrieval, has previously been addressed in various contexts 
in the optics literature, and future work on dynainic spectrum de- 
composition should make full use of that resource. Our particular 
application corresponds to the^problem of retrieving a "complex- 
valued object" (in our case, U) from the modulus of its Fourier 
transforin (|(/| = \/7); this is recognised as a difficult problem 
(McBride, O'Leary and Allen 2004). An acceptable model should 
yield noise-like residuals (in both dynamic and secondary spectra), 
and should satisfy the constraint that the number of free parameters 
be very much less than the number of independent measurements 
of the dynamic spectrum. Our current algorithm fails on the first 
of these criteria. We expect that a globally optimised solution - 
in which all free parameters are simultaneously adjusted - would 
come much closer to achieving noise-like residuals, but we have not 
yet demonstrated this. With such a large number of free parameters, 
a global optimisation would be computationally challenging. 

An important limitation of the algorithm we have presented is 
that it is restricted to finding scattered waves with r ^ 0. Under 
many circumstances this limitation does not cause problems. How- 
ever, if there are multiple refracted images present then problems 
can arise. To see why, we need to consider the properties of the 
starting model: a single component of unit amplitude at r, tj = 0; 
what does this component represent? Because the data themselves 
only carry information about the relative Doppler-shifts, delays and 
phases amongst the various wave components, the choice of ori- 
gin is arbitrary. In practice, because the starting model places all 
of the flux in a single wave, the origin actually corresponds to the 
component which contains the largest flux. If there is only a sin- 
gle refracted image present, then this component will also be the 
component with the smallest delay. However, if there are multiple 
refracted images present, then the brightest of these might well not 
be the path with the smallest delay. In this circumstance the algo- 
rithm will fail to find the scattered and refracted wave components 
which have smaller delays than the starting model - because they 
have r < - and cannot be expected to return an accurate model of 
the electric fields, regardless of how many iterations are performed. 

Our description of the received signal in terms of scattered 
wave components may be used for a variety of purposes; to date 
we have recognised three main applications, which we describe in 
the following sections. 

5.1 Imaging the ionised ISM 

The scattered wave components are identified by the values of their 
delay and Doppler-shift, r and uj. If the scattering occurs in a sin- 
gle, thin screen - as often seems to be the case when organised 
patterns are seen in a dynamic spectrum (Stinebring et al 2001) - 
then there is a direct relationship between these coordinates and 
the apparent positions of the scattering centres. This relationship 
takes a simple form in cases where the observed delays are domi- 
nated by the geometric path delay, leading to parabolic features in 
the secondary spectrum; this circumstance is quite common (Stine- 
bring et al 2001; Cordes et al 2005). If the mapping between (r, u>) 
and position can be determined, then the scattered waves can be 
remapped to give an image of the refracting/scattering centres. 
Such a picture would be complex (containing both ampltitude and 
phase information), and detailed, and should provide valuable in- 
sights into the ionised component of the interstellar medium. In 
particular it will be helpful in elucidating the nature of the anoma- 
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lous scattering and refracting screens wiiich are known to cause 
a variety of phenomena - Extreme Scattering Events (Fiedler et 
al 1987, 1994); refraction and multiple imaging in pulsars (e.g. 
Hewish 1980; Cordes and Wolszczan 1986; Rickett, Lyne and 
Gupta 1994; Hill et al 2005); and intra-day variability in compact 
radio quasars (Kedziora-Chudczer et al 1997; Dennett-Thorpe and 
de Bruyn 2000; BignaU et al 2003). These phenomena are poorly 
understood (Rickett 1991, 2001). 



5.2 Quantifying pulse TOA errors 

If some of the signal arriving at the telescope is delayed, then a 
measurement of the pulse Time-Of- Arrival (TOA) will be affected 
by this delay. As the scattering geometry necessarily changes from 
epoch to epoch, this is a source of systematic errors in precision 
pulsar timing observations. Such studies are important for tests of 
fundamental physics - the equation of state of matter at nuclear 
densities, for example - for low-frequency gravitational wave de- 
tection, and for precision time-keeping (Foster and Backer 1990). 
Eliminating propagation errors in pulse TOAs would therefore be 
valuable. By analysing dynamic spectra with the technique we have 
presented, we obtain detailed information about the various paths 
from source to observer, allowing us to quantify the effects of 
multi-path propagation on the pulse TOAs. In turn this means that 
we should be able to eliminate this contribution to the systematic 
errors in TOAs if high resolution dynamic spectra are recorded in 
tandem with the timing information. 

We can envisage two different schemes for correcting the scat- 
tering errors in pulse TOAs. The simplest scheme would be to com- 
pute the net delay due to all the known paths, and then subtract this 
delay from the corresponding TOA. The second scheme is much 
more ambitious: if the pulse TOA is determined from off-line re- 
duction of baseband data, then we can, in principle, process those 
data in such a way that the electric fields from the various scat- 
tered paths are coherently recombined with the vmscattered signal. 
This generaUses and extends the concept of coherent de-dispersion 
(Hankins and Rickett 1975), whereby the "filtering" imposed by 
propagating the signal through a cold plasma is precisely removed 
by applying the inverse filter. Coherent de-dispersion is now rou- 
tinely used to process baseband data on radio pulsars. By analogy 
with a phase-conjugate mirror, which eliminates wavefront errors 
by exactly reversing light propagation paths, we can term a coher- 
ent recombination of the scattered signal "Virtual Phase Conjuga- 
tion" (VPC). VPC has the potential to completely remove the signal 
filtering imposed by (multi-path) interstellar propagation, and thus 
to eliminate the associated contributions to pulse TOA errors. By 
the same token, VPC should permit studies of pulsar microstruc- 
ture (e.g. Hankins 1996) at very high time-resolution. 



5.3 High resolution imaging of the source 

In §5.1 we noted that the scattered wave decomposition yields, 
fairly directly, an image of the scattering medium. With some fur- 
ther development it should also permit interferometric imaging of 
the pulsar itself. To see why, we need only recall how terrestrial ra- 
dio interferometers operate: the correlation between pairs of signals 
is evaluated as a function of baseline, i.e. separation between the 
antennas. A high visibility amplitude for a given source means that 
the coherent patch is larger than the baseline length, so the source 
must be small. Conversely, if the visibility amplitude is small then 
the baseline is longer than the size of the coherent patch, and we 



have resolved the source on this baseline. Exactly the same con- 
siderations apply to the interference between the scattered waves 
discussed here; in our case the interferometric baseline is simply 
the transverse separation of the paths at the location of the scatter- 
ing screen. 

In the particular approach we have described in this paper the 
scattered waves are identified under the assumption of a point-like 
source (see §2) with an infinitely large coherent patch, so that all 
pairs of scattered waves, j,k yield fringes of amplitude |Sj||Sfc|. 
For very large scattering angles, this approximation must break 
down, and in this case the data can be used to image the source in 
a manner directly analogous to terrestrial radio astronomy, namely 
by quantifying the fringe visibility as a function of baseline length. 
This technique is fundamentally similar to previous investigations 
which used interstellar scattering to constrain the size of pulsar 
radio-emission regions (Gwinn 2001; Gwinn et al 1997; Wolszczan 
and Cordes 1987). There are two main advantages of the method 
proposed here. First, it makes use of a lot more information — all 
the contributing propagation paths are elucidated, thus permitting 
many more constraints on the model brightness distribution of the 
pulsar. Secondly, by knowing the geometry of all the contributing 
paths we can re-order the visibility data into the usual coordinate 
system (the "it, v plane") used for radio astronomical imaging with 
terrestrial interferometers. From that position the many powerful 
concepts, techniques and tools developed for synthesis imaging can 
be brought to bear on the problem of imaging via interstellar scat- 
tering. 

5.4 Interstellar holography 

B.J. Rickett (personal communication, 2003) has previously noted 
the close analogy between his scattered wave solution, derived in 
the weak scattering limit (§2.1), and Gabor holography. The scat- 
tered wave decomposition described in this paper goes beyond the 
weak scattering approximation, but the relationship to holography 
remains strong. In both cases the dynamic spectrum can be re- 
garded as an in-line (Gabor) hologram of the interstellar medium. 
In the weak scattering limit the hologram is recorded with a plane 
reference wave, in effect - the strong, unscattered wave - and can 
be reconstructed with such a wave to yield an image of the scat- 
tering medium. In this paper we have considered the more general 
circumstance where the amplitudes of the scattered components are 
not negligibly small, and their mutual interference must be taken 
into account. In this case the dynamic spectrum is analogous to a 
hologram recorded with an aberrated reference wave, and recon- 
structing an image is no longer so straightforward. Our approach 
will permit images to be reconstructed computationally for this sit- 
uation. The connection with Gabor holography is worth bearing in 
mind because many powerful techniques have been developed in 
that domain, and some may be useful in interstellar holography. 



6 CONCLUSIONS 

We have considered the possibility of representing pulsar dynamic 
spectra in terms of an identifiable collection of electric field com- 
ponents whose mutual interference yields the observed intensity 
structure. Such a representation is of interest because of the di- 
rect link between the properties of the field components and the 
characteristics of the contributing propagation paths through the in- 
terstellar medium. An algorithm for achieving this decomposition 
has been derived, implemented, and tested on high-quaUty data; it 
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works well, although the resulting model is not formally an accept- 
able representation of the data as the residuals are not noise-like. 
Further development of this technique should permit insights into 
a variety of issues relating to interstellar wave propagation and the 
physics of pulsars. 
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